This notebook does some initial exploration of clustering in two
libraries from SCPCP000015: SCPCL000822 and
SCPCL000824.
The clusters are then compared to the results from running
SingleR in the aucell-singler-annotation.sh
workflow. We then look at marker gene expression across the clusters and
assign each cluster to a cell type.
suppressPackageStartupMessages({
# load required packages
library(SingleCellExperiment)
library(ggplot2)
})
# Set default ggplot theme
theme_set(
theme_bw()
)
# set seed
set.seed(2024)
# The base path for the OpenScPCA repository, found by its (hidden) .git directory
repository_base <- rprojroot::find_root(rprojroot::is_git_root)
# The current data directory, found within the repository base directory
data_dir <- file.path(repository_base, "data", "current", "SCPCP000015")
# The path to this module
module_base <- file.path(repository_base, "analyses", "cell-type-ewings")
# path to marker genes file
marker_genes_file <- file.path(module_base, "references", "visser-all-marker-genes.tsv")
# source in helper functions for make_jaccard_matrix() and jaccard()
jaccard_functions <- file.path(module_base, "scripts", "utils", "jaccard-functions.R")
source(jaccard_functions)
# source in helper functions: plot_density() and calculate_sum_markers()
validation_functions <- file.path(module_base, "scripts", "utils", "tumor-validation-helpers.R")
source(validation_functions)
# get louvain, jaccard clusters for a specified value of k (nearest neighbors)
get_clusters <- function(pcs, k) {
clusters <- bluster::clusterRows(
pcs,
bluster::NNGraphParam(
k = k,
type = "jaccard",
cluster.fun = "louvain"
)
)
return(clusters)
}
# define a function to perform clustersweep and get clusters across multiple values of k
cluster_sweep <- function(sce) {
# first perform clustering across parameters
cluster_results <- bluster::clusterSweep(reducedDim(sce, "PCA"),
bluster::NNGraphParam(),
k = as.integer(seq(5, 40, 5)),
cluster.fun = "louvain",
type = "jaccard"
)
# turn results into a data frame
cluster_df <- cluster_results$clusters |>
as.data.frame() |>
# add barcode column
dplyr::mutate(barcodes = colnames(sce)) |>
# combine all cluster results into one column
tidyr::pivot_longer(
cols = ends_with("jaccard"),
names_to = "params",
values_to = "cluster"
) |>
# separate out parameters, nn, function, and type into their own columns
dplyr::mutate(
nn_param = stringr::word(params, 1, sep = "_") |>
stringr::str_replace("k.", "k_"),
cluster_fun = stringr::word(params, 2, sep = "_") |>
stringr::str_remove("cluster.fun."),
cluster_type = stringr::word(params, -1, sep = "_") |>
stringr::str_remove("type.")
) |>
# remove combined params column
dplyr::select(-params)
return(cluster_df)
}
# cluster statistics functions
# get silhouette width and cluster purity for each cluster
# calculates values across all nn_param options used to determine clustering
# all_cluster_results must have nn_param column
get_cluster_stats <- function(sce,
id,
all_cluster_results) {
pcs <- reducedDim(sce, "PCA")
# split clustering results by library id and param used
split_clusters <- all_cluster_results |>
dplyr::filter(library_id == id) |>
split(all_cluster_results$nn_param)
# for each nn_param get cluster width and purity
all_stats_df <- split_clusters |>
purrr::map(\(df){
sil_df <- bluster::approxSilhouette(pcs, df$cluster) |>
as.data.frame() |>
tibble::rownames_to_column("barcodes")
purity_df <- bluster::neighborPurity(pcs, df$cluster) |>
as.data.frame() |>
tibble::rownames_to_column("barcodes")
# join into one data frame to return
stats_df <- sil_df |>
dplyr::left_join(purity_df, by = "barcodes")
return(stats_df)
}) |>
dplyr::bind_rows(.id = "nn_param")
return(all_stats_df)
}
# calculate cluster stability for a single set of clusters using ari
# bootstrap and get ari for clusters compared to sampled clusters
# re-clusters and gets ari across 20 iterations
get_ari <- function(pcs,
clusters,
k) {
ari <- c()
for (iter in 1:20) {
# sample cells with replacement
sample_cells <- sample(nrow(pcs), nrow(pcs), replace = TRUE)
resampled_pca <- pcs[sample_cells, , drop = FALSE]
# perform clustering on sampled cells
resampled_clusters <- get_clusters(resampled_pca, k)
# calculate ARI between new clustering and original clustering
ari[iter] <- pdfCluster::adj.rand.index(resampled_clusters, clusters[sample_cells])
}
ari_df <- data.frame(
ari = ari,
k_value = k
)
}
# get cluster stability for each nn_param cluster results are available for
get_cluster_stability <- function(sce,
id,
all_cluster_results) {
pcs <- reducedDim(sce, "PCA")
# split clustering results by library id and param used
cluster_df_list <- all_cluster_results |>
dplyr::filter(library_id == id) |>
split(all_cluster_results$nn_param)
# for each parameter, get ari values
cluster_stability_df <- cluster_df_list |>
purrr::imap(\(df, k_value){
# make sure k is numeric and remove extra k_
k <- stringr::str_remove(k_value, "k_") |>
as.numeric()
get_ari(pcs, df$cluster, k)
}) |>
dplyr::bind_rows()
return(cluster_stability_df)
}
# plot individual stats for clusters, either purity or width
plot_cluster_stats <- function(all_stats_df,
stat_column,
id) {
plot_df <- all_stats_df |>
dplyr::filter(library_id == id)
ggplot(plot_df, aes(x = nn_param, y = {{ stat_column }})) +
# ggforce::geom_sina(size = .2) +
ggbeeswarm::geom_quasirandom(method = "smiley", size = 0.1) +
stat_summary(
aes(group = nn_param),
color = "red",
# median and quartiles for point range
fun = "median",
fun.min = function(x) {
quantile(x, 0.25)
},
fun.max = function(x) {
quantile(x, 0.75)
}
) +
labs(title = id)
}
# heatmap comparing cluster assignments to SingleR labels
cluster_celltype_heatmap <- function(cluster_classification_df,
id) {
# get a jaccard mtx for each cluster param
jaccard_df_list <- cluster_classification_df |>
dplyr::filter(library_id == id) |>
split(cluster_classification_df$nn_param) |>
purrr::map(\(df) {
make_jaccard_matrix(
df,
"cluster",
"singler_lumped"
)
})
# turn into heatmap list
make_heatmap_list(jaccard_df_list, column_title = glue::glue("{id}-clusters"), legend_match = "k_5", cluster_rows = FALSE)
}
# Density plot looking at marker gene expression across all clusters
# each panel is from a different marker gene list
# each row shows marker gene expression for that cluster
plot_marker_genes <- function(cluster_classification_df,
k_value,
id) {
# pick clustering to use and select those columns
final_clusters_df <- cluster_classification_df |>
dplyr::filter(nn_param == k_value) |>
dplyr::filter(library_id == id)
# grab columns that contain marker gene sums
marker_gene_columns <- colnames(final_clusters_df)[which(endsWith(colnames(final_clusters_df), "_sum"))]
# create individual density plots and combine into one
marker_gene_columns |>
purrr::map(\(column){
plot_density(
final_clusters_df,
column,
"cluster"
)
}) |>
patchwork::wrap_plots(ncol = 1) + patchwork::plot_annotation(glue::glue("{id}-{k_value}-clusters"))
}
# sample and library ids
sample_ids <- c("SCPCS000490", "SCPCS000492")
library_ids <- c("SCPCL000822", "SCPCL000824")
# define input sce files
sce_file_names <- glue::glue("{library_ids}_processed.rds")
sce_files <- file.path(data_dir, sample_ids, sce_file_names) |>
purrr::set_names(library_ids)
# define classification files
singler_results_dir <- file.path(module_base, "results", "aucell_singler_annotation")
classification_file_names <- glue::glue("{library_ids}_singler-classifications.tsv")
classification_files <- file.path(singler_results_dir, sample_ids, classification_file_names) |>
purrr::set_names(library_ids)
# define output files
cluster_results_dir <- file.path(module_base, "results", "clustering", sample_ids)
cluster_results_dir |>
purrr::walk(fs::dir_create)
cluster_file_names <- glue::glue("{library_ids}_cluster-results.tsv")
cluster_output_files <- file.path(cluster_results_dir, cluster_file_names) |>
purrr::set_names(library_ids)
# read in both SCEs
sce_list <- sce_files |>
purrr::map(readr::read_rds)
# read in classification data frame with SingleR results and combine into one
classification_df <- classification_files |>
purrr::map(readr::read_tsv) |>
dplyr::bind_rows(.id = "library_id")
## Rows: 4290 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): barcodes, singler_ontology, singler_annotation, aucell_annotation
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 7414 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (4): barcodes, singler_ontology, singler_annotation, aucell_annotation
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# read in marker genes table
marker_genes_df <- readr::read_tsv(marker_genes_file, show_col_types = FALSE) |>
# account for genes being from multiple sources
dplyr::select(cell_type, ensembl_gene_id, gene_symbol) |>
dplyr::distinct()
# get list of all cell types found
cell_types <- unique(marker_genes_df$cell_type)
# get the sum of expression of all genes for each cell type
gene_exp_df <- sce_list |>
purrr::map(\(sce) {
cell_types |>
purrr::map(\(type){
calculate_sum_markers(marker_genes_df, sce, type)
}) |>
purrr::reduce(dplyr::inner_join, by = "barcodes")
}) |>
dplyr::bind_rows(.id = "library_id")
# join sum expression columns with classification df
classification_df <- classification_df |>
dplyr::left_join(gene_exp_df, by = c("barcodes", "library_id"))
Below we perform Louvain, Jaccard clustering, varying k.
The minimum k is 5 and the maximum k is 40
with a step size of 5.
# perform clustering across k = 5-40 with increments of 5
all_cluster_results <- sce_list |>
purrr::map(cluster_sweep) |>
dplyr::bind_rows(.id = "library_id") |>
dplyr::mutate(
nn_param = forcats::fct_relevel(nn_param, "k_5", after = 0)
)
# get umap embeddings and combine into a data frame with cluster assignments
umap_df <- sce_list |>
purrr::map(\(sce){
df <- sce |>
scuttle::makePerCellDF(use.dimred = "UMAP") |>
# replace UMAP.1 with UMAP1 and get rid of excess columns
dplyr::select(barcodes, UMAP1 = UMAP.1, UMAP2 = UMAP.2)
}) |>
dplyr::bind_rows(.id = "library_id") |>
dplyr::left_join(all_cluster_results, by = c("library_id", "barcodes"))
Below we visualize the cluster assignments using each parameter on a UMAP. This isn’t particularly useful because of clusters, but it can be helpful to see any values of K that may have obvious over or under clustering.
# look at clustering results for each library across all params
library_ids |>
purrr::map(\(id){
umap_df |>
dplyr::filter(library_id == id) |>
ggplot(aes(x = UMAP1, y = UMAP2, color = cluster)) +
geom_point(alpha = 0.5, size = 0.1) +
facet_wrap(vars(nn_param)) +
theme(
aspect.ratio = 1
) +
labs(title = id) +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 1.5)))
})
## [[1]]
##
## [[2]]
Below we calculate a series of statistics:
# get a combined stats dataframe with purity and width for all clusters
all_stats_df <- sce_list |>
purrr::imap(\(sce, library_id){
get_cluster_stats(sce, library_id, all_cluster_results)
}) |>
dplyr::bind_rows(.id = "library_id") |>
dplyr::mutate(
# make sure the order is correct
nn_param = forcats::fct_relevel(nn_param, "k_5", after = 0)
)
# silhouette width for different params
library_ids |>
purrr::map(\(id){
plot_cluster_stats(all_stats_df, width, id)
})
## [[1]]
##
## [[2]]
# cluster purity for different params
library_ids |>
purrr::map(\(id){
plot_cluster_stats(all_stats_df, purity, id)
})
## [[1]]
##
## [[2]]
Here we see that both silhouette width and cluster purity peak around
k=20 or k=25. For SCPCL000822, increasing k after 25 does
not change the width or purity, For SCPCL000824, increasing
k after 20-25 decreases the width and keeps purity around the same.
# calculate cluster stability
stability_df <- sce_list |>
purrr::imap(\(sce, id){
get_cluster_stability(sce, id, all_cluster_results)
}) |>
dplyr::bind_rows(.id = "library_id") |>
dplyr::mutate(
k_value = as.factor(k_value),
# make sure that k = 5 comes first
k_value = forcats::fct_relevel(k_value, "5", after = 0)
)
# plot stability across all values of k
ggplot(stability_df, aes(x = k_value, y = ari)) +
geom_jitter(width = 0.1) +
facet_grid(rows = vars(library_id)) +
labs(title = "Cluster stability") +
stat_summary(
aes(group = k_value),
color = "red",
# median and quartiles for point range
fun = "median",
fun.min = function(x) {
quantile(x, 0.25)
},
fun.max = function(x) {
quantile(x, 0.75)
}
)
Here we see that cluster stability increases as k increases for both libraries. We see a plateau start to appear around k = 20 and k = 25.
SingleRNow we will compare the clustering results to the cell type
assignments obtained from SingleR in the
aucell-singler-annotation.sh workflow. To compare results
we will calculate the Jaccard similarity index between clusters and cell
types. We expect that good clustering will line up with the cell type
assignments so that there is ~ 1 cell type per cluster.
For the plots, we will only display the top 7 cell types and all
other cells will be lumped together into
All remaining cell types.
cluster_classification_df <- all_cluster_results |>
# add in classifications from singler
dplyr::left_join(classification_df, by = c("library_id", "barcodes")) |>
dplyr::mutate(
# first grab anything that is tumor and label it tumor
# NA should be unknown
singler_annotation = dplyr::case_when(
stringr::str_detect(singler_annotation, "tumor") ~ "tumor",
is.na(singler_annotation) ~ "unknown", # make sure to separate out unknown labels
.default = singler_annotation
) |>
forcats::fct_relevel("tumor", after = 0),
# get the top cell types for plotting later
singler_lumped = singler_annotation |>
forcats::fct_lump_n(7, other_level = "All remaining cell types", ties.method = "first") |>
forcats::fct_infreq() |>
forcats::fct_relevel("All remaining cell types", after = Inf),
nn_param = forcats::fct_relevel(nn_param, "k_5", after = 0)
)
# get heatmap for each library showing cluster assignments vs. singler assignments
library_ids |>
purrr::map(\(id){
cluster_celltype_heatmap(
cluster_classification_df,
id
)
})
## Warning in split.default(x = seq_len(nrow(x)), f = f, drop = drop, ...): data length is not a multiple of split variable
## Warning in split.default(x = seq_len(nrow(x)), f = f, drop = drop, ...): data length is not a multiple of split variable
## [[1]]
##
## [[2]]
From these plots, it’s obvious that K = 5 is over clustering and cell types are spread across a lot of clusters. However, for the other values of K, there don’t seem to be any striking differences that are obvious by eye.
We will also look at the expression of marker genes across each cluster. In these plots, each row shows the distribution of the specified marker genes in that cluster. Each panel is labeled with the marker gene expression being plotted.
Based on the above results, we will only show these plots for a value of K = 25.
# plot to look at marker gene expression for all cell types across all clusters
# choose k=25 based on previous stats
library_ids |>
purrr::map(\(id) {
plot_marker_genes(cluster_classification_df,
k_value = "k_25",
id
)
})
## [[1]]
##
## [[2]]
For both samples we see very clear separation between clusters that express tumor marker genes and each of the normal cell marker genes.
Finally, we will assign cell types to all cells in a cluster based on a cluster having high expression of a set of marker genes.
SCPCL000822 that has high
mesenchymal gene expression and is predominantly chondrocytes so we will
label this cell as “chondrocytes”.SCPCL000824.# get umap embeddings for plotting later
umap_df <- umap_df |>
dplyr::select(library_id, barcodes, UMAP1, UMAP2) |>
dplyr::distinct()
# assign clusters based on marker gene expression and cell types that are prominent in that cluster
# use the jaccard matrices to help guide assignments
cluster_umap_df <- cluster_classification_df |>
dplyr::left_join(umap_df, by = c("library_id", "barcodes")) |>
dplyr::filter(nn_param == "k_25") |>
dplyr::mutate(
cluster_celltype = dplyr::case_when(
library_id == "SCPCL000822" & cluster %in% c(1, 2, 4, 5, 9) ~ "tumor",
library_id == "SCPCL000822" & cluster == 7 ~ "endothelial cell",
# cluster 3 has predominantly chondrocytes and has MSC marker gene expression
library_id == "SCPCL000822" & cluster == 3 ~ "chondrocyte",
# clusters with high "immune" marker gene expression are macrophage
library_id == "SCPCL000822" & cluster == 6 ~ "macrophage",
# cluster 8 has high MSC marker gene expression and a mix of MSC like cell types
library_id == "SCPCL000822" & cluster == 8 ~ "mesenchymal-like cell",
library_id == "SCPCL000824" & cluster %in% c(1, 2, 3, 4, 5, 6, 11) ~ "tumor",
library_id == "SCPCL000824" & cluster == 8 ~ "endothelial cell",
# cluster 10 has high immune markers and is mostly macrohpages
library_id == "SCPCL000824" & cluster == 10 ~ "macrophage",
# cluster 7 has high MSC marker gene expression and a mix of MSC like cell types
library_id == "SCPCL000824" & cluster == 7 ~ "mesenchymal-like cell",
# if they didn't have high expression of any markers, we label them as unknown
.default = "unknown"
),
# make tumor first
cluster_celltype = forcats::fct_relevel(cluster_celltype, "tumor", after = 0)
)
# look at cell type assignments on the UMAP
ggplot(cluster_umap_df, aes(x = UMAP1, y = UMAP2, color = cluster_celltype)) +
geom_point(alpha = 0.5, size = 0.5) +
facet_grid(cols = vars(library_id)) +
theme(
aspect.ratio = 1,
legend.position = "bottom"
) +
guides(color = guide_legend(override.aes = list(alpha = 1, size = 1.5))) +
scale_color_brewer(palette = "Dark2")
We can also validate that these cells do in fact have higher expression of the marker genes as compared to other cells by plotting the marker gene expression on UMAPs.
# grab columns that contain marker gene sums
marker_gene_columns <- colnames(cluster_umap_df)[which(endsWith(colnames(cluster_umap_df), "_sum"))]
# for each sample show a faceted umap by cell type
# color is gene expression of all marker genes for that cell type
library_ids |>
purrr::map(\(id){
plot_df <- cluster_umap_df |>
dplyr::filter(library_id == id) |>
tidyr::pivot_longer(
cols = all_of(marker_gene_columns),
names_to = "cell_type",
values_to = "gene_expression"
) |>
dplyr::mutate(cell_type = stringr::str_remove(cell_type, "_sum"))
ggplot(plot_df, aes(x = UMAP1, y = UMAP2, color = gene_expression)) +
geom_point(alpha = 0.5, size = 0.2) +
facet_wrap(vars(cell_type)) +
theme(aspect.ratio = 1) +
scale_color_viridis_c() +
labs(
title = id
)
})
## [[1]]
##
## [[2]]
Here we will save the assigned clusters and cell type annotations for each library. Each file will have the following columns:
barcodes: Unique cell barcodecluster: Assigned clustercluster_celltype: Cell type assigned based on
clustering and marker gene expressionnn_param: Number of nearest neighbors used as a
parameter for clusteringsingler_ontology: Cell ontology identifier for assigned
cell type by SingleRsingler_annotation: Cell type assigned by
SingleRconsensus_celltype: Consensus between cluster
assignments and SingleR assignments for tumor cells. Here
any cell in a cluster with tumor cells will be labeled with
tumor and all other cells will be labeled with the
SingleR annotation.cluster_output_files |>
purrr::iwalk(\(file, id){
cluster_umap_df |>
dplyr::filter(library_id == id) |>
dplyr::select(barcodes, cluster, cluster_celltype, nn_param, singler_ontology, singler_annotation) |>
dplyr::mutate(consensus_celltype = dplyr::if_else(cluster_celltype == "tumor", "tumor", singler_annotation)) |>
readr::write_tsv(file)
})
# record the versions of the packages used in this analysis and other environment information
sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.6.7
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] ggplot2_3.5.1 SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0 Biobase_2.64.0 GenomicRanges_1.56.1
## [6] GenomeInfoDb_1.40.1 IRanges_2.38.1 S4Vectors_0.42.1 BiocGenerics_0.50.0 MatrixGenerics_1.16.0
## [11] matrixStats_1.3.0
##
## loaded via a namespace (and not attached):
## [1] rlang_1.1.4 magrittr_2.0.3 clue_0.3-65 GetoptLong_1.0.5 pdfCluster_1.0-4
## [6] compiler_4.4.0 DelayedMatrixStats_1.26.0 png_0.1-8 vctrs_0.6.5 stringr_1.5.1
## [11] shape_1.4.6.1 pkgconfig_2.0.3 crayon_1.5.3 fastmap_1.2.0 XVector_0.44.0
## [16] labeling_0.4.3 scuttle_1.14.0 magic_1.6-1 utf8_1.2.4 rmarkdown_2.27
## [21] tzdb_0.4.0 UCSC.utils_1.0.0 ggbeeswarm_0.7.2 bit_4.0.5 purrr_1.0.2
## [26] xfun_0.46 bluster_1.14.0 cachem_1.1.0 zlibbioc_1.50.0 beachmat_2.20.0
## [31] jsonlite_1.8.8 highr_0.11 DelayedArray_0.30.1 BiocParallel_1.38.0 tweenr_2.0.3
## [36] parallel_4.4.0 cluster_2.1.6 R6_2.5.1 bslib_0.8.0 RColorBrewer_1.1-3
## [41] stringi_1.8.4 jquerylib_0.1.4 Rcpp_1.0.13 iterators_1.0.14 knitr_1.48
## [46] readr_2.1.5 Matrix_1.7-0 igraph_2.0.3 tidyselect_1.2.1 rstudioapi_0.16.0
## [51] abind_1.4-5 yaml_2.3.10 doParallel_1.0.17 codetools_0.2-20 lattice_0.22-6
## [56] tibble_3.2.1 withr_3.0.0 evaluate_0.24.0 polyclip_1.10-7 circlize_0.4.16
## [61] pillar_1.9.0 BiocManager_1.30.23 renv_1.0.7 foreach_1.5.2 geometry_0.4.7
## [66] generics_0.1.3 vroom_1.6.5 rprojroot_2.0.4 hms_1.1.3 sparseMatrixStats_1.16.0
## [71] munsell_0.5.1 scales_1.3.0 glue_1.7.0 tools_4.4.0 BiocNeighbors_1.22.0
## [76] forcats_1.0.0 fs_1.6.4 Cairo_1.6-2 grid_4.4.0 tidyr_1.3.1
## [81] colorspace_2.1-1 GenomeInfoDbData_1.2.12 patchwork_1.2.0 beeswarm_0.4.0 ggforce_0.4.2
## [86] vipor_0.4.7 cli_3.6.3 fansi_1.0.6 viridisLite_0.4.2 S4Arrays_1.4.1
## [91] ComplexHeatmap_2.20.0 dplyr_1.1.4 gtable_0.3.5 sass_0.4.9 digest_0.6.36
## [96] SparseArray_1.4.8 rjson_0.2.21 farver_2.1.2 htmltools_0.5.8.1 lifecycle_1.0.4
## [101] httr_1.4.7 GlobalOptions_0.1.2 bit64_4.0.5 MASS_7.3-61